
clc; clear; close all;

N = 500000;  % 迭代步数

w_QIM = 0.8; 
w_SIM = 0.5;
w_TGM = 0.8;
w_CGM = 0.2;

% 设置每个模型的幅值数组 A
A_array_QIM = [1.0, 1.2, 1.4, 1.6];
A_array_SIM = [0.1, 0.2, 0.3, 0.4];
A_array_TGM = [1.0, 1.2, 1.4, 1.6];
A_array_CGM = [1.0, 1.2, 1.4, 1.6];

% 颜色列表 (可根据喜好增减)
color_list = {'r','b','g','k','m','c'};

%% (a) Q-IM
subplot(2,2,1);
hold on; grid on; box on;
legend_str = cell(length(A_array_QIM),1);
for idx = 1:length(A_array_QIM)
    A_val = A_array_QIM(idx);
    [v, i] = QIM(w_QIM, N, A_val);
    plot(v, i, color_list{idx}, 'LineWidth',1.2);  % 横轴 i，纵轴 v
    legend_str{idx} = sprintf('A=%.1f', A_val);
end
xlabel('v_1(n), A'); ylabel('i_1(n), V');
title('(a) Q-IM');
legend(legend_str, 'Location','best','NumColumns',2);

%% (b) S-IM
subplot(2,2,2);
hold on; grid on; box on;
legend_str = cell(length(A_array_SIM),1);
for idx = 1:length(A_array_SIM)
    A_val = A_array_SIM(idx);
    [v, i] = SIM(w_SIM, N, A_val);
    plot(v, i, color_list{idx}, 'LineWidth',1.2);
    legend_str{idx} = sprintf('A=%.1f', A_val);
end
xlabel('v_2(n), V'); ylabel('i_2(n), A');
title('(b) S-IM');
legend(legend_str, 'Location','best','NumColumns',2);

%% (c) T-GM
subplot(2,2,3);
hold on; grid on; box on;
legend_str = cell(length(A_array_TGM),1);
for idx = 1:length(A_array_TGM)
    A_val = A_array_TGM(idx);
    [v, i] = TGM(w_TGM, N, A_val);
    plot(v, i, color_list{idx}, 'LineWidth',1.2);
    legend_str{idx} = sprintf('A=%.1f', A_val);
end
xlabel('v_3(n), V'); ylabel('i_3(n), A');
title('(c) T-GM');
legend(legend_str, 'Location','best','NumColumns',2);

%% (d) C-GM
subplot(2,2,4);
hold on; grid on; box on;
legend_str = cell(length(A_array_CGM),1);
for idx = 1:length(A_array_CGM)
    A_val = A_array_CGM(idx);
    [v, i] = CGM(w_CGM, N, A_val);
    plot(v, i, color_list{idx}, 'LineWidth',1.2);
    legend_str{idx} = sprintf('A=%.1f', A_val);
end
xlabel('v_4(n), V'); ylabel('i_4(n), A');
title('(d) C-GM');
legend(legend_str, 'Location','best','NumColumns',2);

sgtitle('四种离散忆阻器在不同幅值下的钳制回滞曲线');

%% 封装函数：Q-IM 模型
function [v_array, i_array] = QIM(omega, N, A)
    % 初始状态
    phi = 0;
    v_array = zeros(1, N);
    i_array = zeros(1, N);
    for n = 1:N
        v = A * sin(omega * n);
        i = (phi^2 - 2) * v;
        v_array(n) = v;
        i_array(n) = i;
        % 更新 phi
        phi = phi + v;
    end
end

%% 封装函数：S-IM 模型
function [v_array, i_array] = SIM(omega, N, A)
    phi = 0;
    v_array = zeros(1, N);
    i_array = zeros(1, N);
    for n = 1:N
        v = A * sin(omega * n);
        i = sin(pi * phi) * v;
        v_array(n) = v;
        i_array(n) = i;
        phi = phi + v;
    end
end

%% 封装函数：T-GM 模型
function [v_array, i_array] = TGM(omega, N, A)
    phi = 0;
    v_array = zeros(1, N);
    i_array = zeros(1, N);
    for n = 1:N
        v = A * sin(omega * n);
        i = sin(phi) * v;
        v_array(n) = v;
        i_array(n) = i;
        % phi 更新：带有非线性项
        phi =  0.12 * phi * abs(phi) + 0.86 * phi - 0.02 * (phi^3) + 0.1 * v;
    end
end

%% 封装函数：C-GM 模型
function [v_array, i_array] = CGM(omega, N, A)
    phi = 0;
    v_array = zeros(1, N);
    i_array = zeros(1, N);
    for n = 1:N
        v = A * sin(omega * n);
        i = cos(phi) * v;
        v_array(n) = v;
        i_array(n) = i;
        % phi 更新：带 tanh 非线性项
        phi = 0.1 * tanh(phi^3) + 0.99 * phi + 0.1 * v;
    end
end
